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Abstract 

A Monte Carlo simulation of chemotactic bacteria is developed on the basis 
of the kinetic model and is applied to a one-dimensional traveling population 
wave in a microchanneh In this simulation, the Monte Carlo method, which 
calculates the run-and-tumble motions of bacteria, is coupled with a hnite 
volume method to calculate the macroscopic transport of the chemical cues 
in the environment. The simulation method can successfully reproduce the 
traveling population wave of bacteria that was observed experimentally and 
reveal the microscopic dynamics of bacterium coupled with the macroscopic 
transports of the chemical cues and bacteria population density. The results 
obtained by the Monte Carlo method are also compared with the asymptotic 
solution derived from the kinetic chemotaxis equation in the continuum limit, 
where the Knudsen number, which is dehned by the ratio of the mean free 
path of bacterium to the characteristic length of the system, vanishes. The 
validity of the Monte Carlo method in the asymptotic behaviors for small 
Knudsen numbers is numerically verihed. 
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1. Introduction 


Due to innovative developments in biotechnology, including tissue engi¬ 
neering and cell engineering, studies on active fluids composed of biological 
entities have recently drawn increasing attention from physical and mathe¬ 
matical scientists. The pattern formations and fluid flows that are sponta¬ 
neously generated by the internal motions of active entities, which interact 
with the local environment, are studied from a physical point of view at both 
the microscopic and macroscopic levels. A suspension of chemotactic bacte¬ 
ria, e.g., E. Coli, is a typical example of the active fluid, in which the bacteria 
create the collective motions, and the macroscopic patterns are created spon¬ 
taneously via interactions with chemical cues, whose local concentrations also 

vary.ll|,y,y,l4| 

The macroscopic transport phenomena created by the collective motion 
of chemotactic bacteria can be modeled using the coupled reaction-diffusion 
equations for nutrients (which are consumed by bacteria), chemoattractants 
(which are secreted by bacteria), and bacterial density. The basic idea of 
utilizing the react ion-diffusion equation to describe the collective motion of 
chemotactic bacteria was first introduced by Keller and Segelp, l6|, and both 
mathematical and physical studies on the Keller-Segel model have accumu¬ 
lated. Although the microscopic basis of the Keller-Segel model was more 
recently laid by the asymptotic analysis of the kinetic chemotaxis model, 
this type of modeling is originally based on a phenomenological point of 
view. The macroscopic transport phenomena are consequences of the mi- 
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grations of bacteria, molecular diffusion of chemical cues (i.e ., nutrients and 
chemoattractants), and their mutual interactions. 8|, l9|, [lO, [ill] Thus, the 
transport phenomena of chemotactic bacteria set an important multiscale 
problem to be resolved from physical and mathematical points of view. The 
mesoscopic model to describe the connection between the microscopic dy¬ 
namics of bacterium and the macroscopic transports of chemical cues and 
bacteria population density takes on major signihcance. 

The kinetic approach to describin g th e collective motion of chemotactic 


bacteria was hrst introduced by Alt 


12| and then further developed jl3 |. 


In the collective motions of the bacteria, each bacterium repeats a simple 
run-and-tumble motion, where they run ballistically in some duration and 
subsequently change their directions randomly. In kinetic modeling, the run- 
and-tumble motion of bacterium is assumed to be a stochastic process, and 
the time evolution of the density of bacteria with a velocity of n, f{t, x, v) 
is described by a variant of the Boltzmann transport equation for gases (al¬ 
though the present kinetic transport equation does not involve the quadratic 
operator for collisions). IJ, [l^ The transition (or scattering) kernel in the 
kinetic chemotaxis equation involves a model response function that stochas¬ 
tically determines the run durations of each bacterium according to the en¬ 
vironment experienced by the bacteria along their pathways. The response 
function takes a major role to reproduce the chemotactic motions of each 
bacterium and their collective behaviors. Thus, the kinetic approach is a 
promising candidate to investigate the connection between the microscopic 
dynamics of bacteria and macroscopic transport phenomena. 

Recently, the studies on the kinetic chemotaxis have been actively per- 
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formed by mathematical scientists. Several researchers have succeeded in 
deriving the macroscopic continuum equations for chemotactic bacteria from 
the kinetic transport equation, and the connection between the macroscopic 
continuum description and the mesoscopic kinetic description has gradually 


been revealed. 


171, [iSl y, |20l l2l|, I 22 , 


23 


m 


25 


26| A comprehensive 


mathematical study on the traveling wave of the kinetic chemotaxis model 


is also given by Galvez. ^ The kinetic chemotaxis model has also been uti¬ 
lized in the analysis of the experimental results. Saragosti et al described the 
largely biased motions of bacteria toward regions with a high concentration 
of chemical ci^, and they proposed a kinetic model based on experimental 
observations. 


The numerical simulations of the kinetic chemotaxis model also takes 
on a signihcance in understanding the multiscale mechanism in the collective 
motions of chemotactic bacteria or analyzing the problems occurring in prac¬ 
tical engineering and biological systems. In the simulation methodologies for 
solving the kinetic chemotaxis models, the difficulty arises in the treatment 
of the response function in the scattering kernel. It should be noted that the 
response function depends not only on the instantaneous spatial distribution 
of chemical cues but on the temporal variations along the running directions 
of each bacterium. Recently, a Cartesian- mesh-based numerical method to 
accurately solve the kinetic chemotaxis equation was developed by Yang and 
Filbet and applied to various one- and two-dimensional problems. 29j In the 
method, an elaborate numerical algorithm is employed to treat the scattering 
kernel of the kinetic transport equation. 

In the present study, a Monte Carlo simulation method for chemotactic 
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bacteria in a three-dimensional system is newly developed on the basis of the 


kinetic chemotaxis equation used in Ref. 


, and the Monte Carlo method 


is applied to the traveling population wave in a microchannel. 


30 


31 


m, 


33| 


Because the Monte Carlo simulation employs a particle-based method, the 
treatment of the response function, which may depend on the memory of 
bacterium along its pathway}^. 1^. 1^ . 1^. 1^ . is simplihed. Some nu¬ 

merical studies toward this direction were also put forward by Rousset and 


Samaey. The microscopic dynamics of bacterium and its relation to 


the macroscopic transports of chemical cues are investigated in detail. The 
simulation results are also compared to the asymptotic solution obtained from 
the kinetic chemotaxis model in order to validate the accuracy of the Monte 
Carlo method in the asymptotic behaviors. Thus, the purpose of the present 
study is to propose a new Monte Carlo simulation method and to demon¬ 
strate the validity and utility of the method numerically. Incidentall y, o ther 
Monte Carlo simulations, e.g., the Brownian dynamics simulations j34l. l35l | 
and the velocity-jump simulation involving a memory kernel for the ki¬ 
netic chemotaxis are also proposed. In the present Monte Carlo method, the 
velocity jump process described by the kinetic chemotaxis equation involving 
a response function of the tumbling rate and a persistence of the reorientation 
angle is utilized. 

In Sec. II, the basic model, i.e., the kinetic chemotaxis model, and its 
non-dimensional form are described. In Sec. Ill, the methodology of the 
present Monte Carlo simulation is explained in detail. In Sec. IV, the nu¬ 
merical results of the Monte Carlo simulations are presented for the macro¬ 
scopic transports, the microscopic dynamics, and the effects of variations in 
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the parameters of the response function. Special focus is placed on the re¬ 
lation between the microscopic dynamics of bacterium and the macroscopic 
transports of the chemical cues and bacteria population density. In Sec. V, 
the results of Monte Carlo simulations are compared to the asymptotic solu¬ 
tion of the kinetic chemotaxis model in the continuum limit to examine the 
accuracy of the Monte Carlo method in the asymptotic behaviors. Finally, 
concluding remarks and an outlook are given in Sec. VI. 


2. Basic Model 


The basic model is described in the reference 


28l |. In the modeling, 


the macroscopic transport of chemical cues, i.e., the nutrients consumed by 
bacteria and the chemoattractants secreted by bacteria, is described by con¬ 
tinuum react ion-diffusion equations, whereas the run-and-tumble motions of 
bacteria are described by a kinetic equation. Because a small bacterium, e.g., 
E. coli, is assumed, the cells are not able to choose directly the preferential 
direction of motion toward the region of a high concentration of chemical 
cues by measuring the head-to-tail gradient of the chemical cues. Instead, 
the cells detect the preferential direction by sensing the temporal variation 
of chemical cues experienced along their pathways. This sensing strategy of 
bacteria is accounted for in the response function in the kinetic equation. 


2.1. Basic equations 

In the following, the basic kinetic equations are explained briefly. f{t, x, v) 
represents the density of bacteria with a velocity n at a time t and a position 
X, and N{t, x) and S{t, x) represent the concentration helds of nutrient and 
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chemoattractant, respectively, at a time t. These quantities are described by 
the following equations: 


and 


dS ^ d^S r, . f rr n, , 


dN f 


( 1 ) 

( 2 ) 


dt dXa 


r(u, v')f{t, X, v')d v' 


Jv'&V 

+ rf{t,x,v). 


T{v\ v)f{t, X, v)d v' 


Iv'GV 


( 3 ) 


Here, Ds and Dat are the diffusion coefficients for the chemoattractant and 
nutrient, respectively, a is the degradation rate of chemoattractant, b is the 
production rate of chemoattractant by a bacterium, and cN is the consump¬ 
tion rate of nutrient by a bacterium. Hereafter, the subscripts a, [3, and 
7 represents the index in Cartesian coordinates, i.e., {q:,/9,7} = {x,y,z}. 
Equation ([3]) is a variant of the Boltzmann equation for gases; in Eq. (|3]), the 
transition (scattering) kernel T{v, v') stands for the tumbling event of a bac¬ 
terium; during this event, the bacterium changes from velocity v' to v. The 
velocity space V is bounded and symmetric. In this study, we solely consider 
the case for bacteria with a preferential velocity Vq. Thus, the integral domain 
V represents the surface of sphere of a radius Vq, i.e., V = = Vq}. 

The last term in Eq. ([3]) represents cell division, where r is the division rate 
of a bacterium (r = In 2 /r 2 , where T 2 is the mean doubling time). 

The transition kernel T{v,v') is assumed to be split into two contribu¬ 
tions; one is the tumbling rate A(ij'), and the other is the reorientation effect 
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during tumbles K{v,v'). 


T{v,v') = X(v')K(v,v'), 


(4) 


with the condition 


K{v, v')dv = 1. 


(5) 


'vGV 


For the tumbling rate A(n), we assume that the bacteria are sensitive to 
the temporal variations of chemical cues along their pathways via logarithmic 
sensing mechanics j^, H] and that the contributions of each chemical cue 
are independent and additive. Then, the tumbling rate A(n) can be written 
as 


HV) = 1 (A„(^') + As(t.')) 

(6a) 

where the is the material derivative along the trajectory with a velocity 
v'. The response functions 'ipN and 'ips are both positive and decreasing 
because bacteria are less likely to tumble (thus perform longer runs) when 
the chemical cues increase. In the present model, the following analytic 
function is chosen as the response function -0; 




D logiV 


Dt 


+ ^s 


D\ogS 


Dt 


'ipiX) = 'ipo — ytanh 



(7) 


where i/jq is the basal mean tumbling frequency, y is the modulation ampli¬ 
tude, and 5“^ is the characteristic time, which represents the stiffness of the 
response function. 

K{v,v') accounts for persistence in the successive run trajectories after 
the tumbles. Q When we consider the uniform scattering kernel, K is & 












constant, i.e., K = and T{v,v') is proportional to the tnmbling rate 

A(n'). Persistence in the snccessive rnns, v' —)■ n, can be described nsing 
the reorientation angle 6, i.e., 


as 


cos 9 = 


K{v, v') oc G 


V ■ V 
^0 


1 — cos 9 




( 8 ) 


(9) 


where a is the standard deviation of reorientation angle 9, i.e., a = \/< 6*^ >, 
and G(X) is an increasing fnnction, and its proportionality constant is de¬ 
termined nsing the normalization condition Eq. ([5]). For G{X) = exp(X), 
one can write the Von Mises distribntion, 

1—cos©'' 

( 10 ) 


exp 


(-i^) 


27rV?(T^ (1 — e 


K{v, v') = 


Experiments also demonstrated that the tnmbling freqnency is strongly cor¬ 
related with the reorientation angle. 


Thns, one can assnme a linear cor¬ 


relation between the standard deviation of reorientation angle a and the 
tnmbling freqnency \{y') as 


a = ai + a2\{v'). 


( 11 ) 


Note that the rotational diffnsivity and long-term memory of the bacteria 
are not involved in the present kinetic chemotaxis model. 


2.2. Non-dimensionalization 

We introdnce the non-dimensional time t, space x, and velocity e, which 
are defined as 

X = x/Lo, e = v/Vo, i = t/to {to = Lq/Vo) . (12) 
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We introduce the characteristic length Lq, which may affect the geometry of 
problem. The bacterial density f{t,x,v) is also non-dimensionalized as 

f{i,x,e) = f{t,x,v)/{po/V^). (13) 


Using the non-dimensional quantities, Eqs. ([I]) - ([2]) are written as fol¬ 
lows. 

rl.Q . ^ ^ r 

f{t,x,e')dn{e'), (14) 


dS d^S 


all e' 


dN - f 

= Dn^-cN f{t,x,e')dQ{e'), 


dt 


dx‘i, 


(15) 


'all e' 


where Dn,s = Dn,s/{L l/to), d = toa, b = 1, and c = potoc. The concentra¬ 
tion of chemoattractant S is scaled by the reference quantity Sq = potob and 
that of the nutrient N is scaled by an arbitrary reference quantity Nq. Note 
that the effect of the initial average density po of bacteria is involved only 
via the non-dimensional parameter c. 

The kinetic chemotaxis equation of bacterial density f{t, x, v) is written 
in non-dimensional form as 


^ + e [ 

Ot 0X(^ Jall e' 


X{e')K{e,e')f{i,x,e')dn{e') - A(e)/(f, :r, e) +f/(f,;E,e), 

(16) 


where A = (Lo/Vo)A, K = V^K, and f = (Lo/Vo)r. We also write A(e') as 
A(e') = '00^(60; where ipo = (Lo/Vo)'0oi and the modulation of the tumble 
frequency 4/(e') is written as 


*(e') = 2 




D log N 


Dt 


+ 


D\ogS 


Dt 


with 


X 


i’s,N(X) = 1 - xs.jvtanh I y I ■ 


(17) 


(18) 
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Here x = x/'^o and <5 = {Lq/Vq)5. Note that ! Lq) corresponds 

to the Knudsen number in the rarehed gas dynamics because the product of 
the bacterial velocity Vq and the inverse of the mean tumble frequency ■, 
which is the mean free time, represents the mean free path. K{e,e') is 
written as 


iX(e,e') = 


exp 




27rcT^ ( 1 


2 


(19) 


K{e,e') also satishes K(e, e')dQ(e) = 1. 

The standard deviation of reorientation angle a can be written as 


a = cTi + a2^(e'), 


( 20 ) 


where the constant values of ai and <72 are determined from the maximum 
and minimum values of the standard deviation of reorientation angle, ciMax 
and (Tmin, as 


(^Max — CTi + Cr2TMax, 

<^min T Cr2^min) 


(21a) 

(21b) 


where ^Max and Tmin are the maximum and minimum values of Eq. fll7p . 
respectively, and are set to Tnax = 1-4 and Tmin = 0.6 for Eqs. flTT)) and 
03 ; these values are estimated from experimental results [^ . 


3. Simulation Method 

In the present paper, we consider the one-dimensional macroscopic trans¬ 
ports of chemical cues, i.e., dS/dy = dS/dz = 0 and dN/dy = dN/dz = 0, 
while the three-dimensional motions of bacteria are calculated by the Monte 
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Carlo method as described below. The spatial domain is discretized into a 
uniform lattice mesh system with ly x cubes of a constant side length 
Ah. Here, Jq, is the number of mesh intervals in the a direction of Cartesian 
coordinates. The macroscopic transports of chemical cues are calculated on 
this lattice mesh system. Because we consider the one-dimensional trans¬ 
ports of chemical cues along x coordinate, we can set ly = Iz = 1. The 
lattice mesh nodes are expressed as Xi {iAx) {i = 0,1, • ■ • ,Ix) and the cen¬ 
ters of each lattice are expressed as (i = 0,1, • • • , — 1). See also Fig. 
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Figure 1: The lattice mesh system and the simulation particles in the lattice site. The 
run-and-tumble motions of particles are calculated by the Monte Carlo method while 
the concentrations of chemical cues are calculated by the finite volume method. The 
specular reflection condition for the particles and the non-flux condition for chemical cues 
are considered at x = xq and x/^ while the periodic conditions are considered for both 
particles and chemical cues in y and z directions. 


Equations flTT|) and flT^ are calculated with using a hnite volume scheme. 



At 


Ax ydx J 



as”+‘ + bill 


n+1 


^ I _ I 

At 


Dn 

Ax 




cvr+vr 


n-l-l 


( 22 ) 

(23) 
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for i = 0,1, • • • , /a; — 1, where the superscript n represents the time step 
number, the subscript i represents the zth lattice site, i.e., F" = F{nAt, £j+i) 
{F = S or N), and {dF/dx)'^ is dehned as 


\dx J ■ Ax 

for i = 1, • • • ,1^ — 1. {dF/dx)o and {dF/dx)i^ are determined by the bound¬ 
ary conditions. The macroscopic density of bacteria in the ith lattice site pi 
is dehned as 



P\ = 


Ax^ 


f{nAi, X, e')dVL{e')d 


' 2 th site J all e' 


(25) 


and is calculated from the number of simulation particles used in the Monte 
Carlo simulation as described in the next paragraph. Note that the mesh 
interval Ax must be sufficiently small to resolve the macroscopic transports 
of chemical cues and the time-step size At must be smaller than the diffusion 
time, i.e.. At < QAAxP'/D s^n- 

The three-dimensional motions of bacteria are calculated in each lattice 
site by the Monte Carlo method. We use a constant and uniform weight wq 
for a single simulation particle; this weight represents the number of bacteria 
corresponding to each simulation particle. That is, the number of bacteria 
in the ith lattice site is written as 


(LoAxfpoPi = wopi, (26) 

where pi is the number of simulation particles contained in the ith lattice 
site. Thus, the macroscopic density at the ffih lattice site pi is calculated from 
the number of simulation particles pi via Eq. fl26|) . Note that the initial total 
number of simulation particles Mq and the weight Wq are not related to the 
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physical results but to the accuracy of the Monte Carlo simulation. The effect 
of the initial mean density of bacteria po is involved via the non-dimensional 
parameter c. 

The Monte Carlo simulation is conducted using the following steps. Here¬ 
after, the position and velocity of the Ith particle are expressed as r(p and 
e(p, respectively. 

0. Att = 0, the simulation particles are distributed according to the initial 
density /°(e) {i = 0, ■■■, — 1). The number of simulation particles 

in the ith lattice site, is determined by Eq. fl2^ . In each lattice site, 
simulation particles are distributed uniformly at random positions and 
their velocities e are determined by the probability density /°(e)/p°. 

1. Particles move with their velocities for a duration At: 

r5« = r5, + e5,A( (i = 1, ■ • • , M"), (27) 

where M” is the total number of simulation particles at time step n. 
The particles that move beyond the boundaries are removed, and new 
ones are inserted according to the boundary conditions. 

The number of particles at each lattice site, (z = 0, • • ■ , Ja: — 1) is 
also counted. 

2. At each lattice site, the concentrations of chemical cues and 

(z = 0, • • • , Ja; — 1) are calculated by Eqs. fl2^ and fl2^ with Eq. fl26|) . 

3. The tumbling of each particle is calculated using the scattering kernel in 
Eq. ffT6|) . The tumbling of the Zth particle may occur with a probability 
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is calculated from the temporal variations of chemical cues along 
the pathway of the Zth particle, —)■ in the successive time steps. 

We write the chemical cues experienced by the Ith particle at the time 
step n as 5^^ and NJ^y Then, 

^ (^^5(0+V'v(/)) , (28) 


with 


i’S{i) = 1 - Xs tanh 
'ipNn) = 1 - XArtanh 


/log^"+'-log^, 


I 


SAi 

j-n+l 


logiV(^+^-logiV, 

SAi 


n 

(0 


(29a) 

(29b) 


and are the local concentration of chemical cues at the position 
of the Ith particle, r^y and are calculated by the interpolation using 
the concentrations of the chemical cues at the lattice site in which the 
particle in involved and the linear gradients between the neighboring 
lattice sites. Thus, when the Ith particle is involved in the ith lattice 
site, {F = S or N) is calculated as 


/ 3F\ ” 


(30) 


where the linear gradient is calculated using the central difference be¬ 
tween neighboring lattice sites as 


dF^ ” 
dx 


i+i 


fdFY fdF\ 


i+i. 


(31) 


\dx J 

where {dF/dx)i is dehned in Eq. fl2T)) . Using the linear interpolation 
Eq. (1501) . a simulation particle that stays at the same lattice site after 
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a single time step passes can sense the chemical gradients along its 
pathway. 

For the particle that is judged to tumble, say the I'th. particle, a new run 
direction after the tumbling, is determined using the probability 

in Eq. ([19]). By using the stochastic variables 6 and 0 

calculated as 


cos 6* 


1 + cr^ log 


e“^ + (1 - e~^)U^ , 


0 = 271112 , 


(32a) 

(32b) 


where Ui and U 2 are uniform random numbers, the new run direction 
of the /'th particle is obtained by 


cos 6 + sin 6 cos (j) + sin 9 sin cj), 


(33) 


where e” and are the normal vectors to and are mutually per¬ 
pendicular. For example, they can be written as 


e 


n 

1 



(34a) 

(34b) 


4. For all simulation particles, divisions are calculated with a uniform 
probability fAt. For a particle that is judged to undergo division, e.g., 
the Ith particle, a new particle with the run direction e(/) is created at 
a random position in the same lattice site. Note that the total number 
of simulation particles is updated as 
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5. Return to step 1 with the obtained r(p, e(/), S^i), and iV(p (/=!,• • • ,M) 
at the new time step. 

4. Results 

The simulation method described in the previous section is applied to 
the traveling population wave of bacteria in a microchanneh The parameter 
values used in the simulations, unless otherwise stated, are summarized in 
Tables [U and [2l 


Table 1: Reference quantities. 
Lq 1 [mm] 

Vb 25 [/tm/s] 

to 40 [s] 

po 5x10® [cell/mL] 


The model parameters are chosen to reproduce the experimental and 


numerical results obtained in Ref. 


28| . The numerical parameters ui and 


(J 2 are determined from Eq. fl2l]) with UMax = 1.5 and (Tmin = 1.3, which also 
correspond to the experimentally reported values ^ 


The total number of simulation particles is ten times larger than the 
number of bacteria in the experimental system, i.e., the weight r(;o=0.1, to 
obtain an accurate solution to Eqs. flTTll - ffTOll by reducing the noise that 
arises in the Monte Carlo method. 

The one-dimensional lattice mesh system is used for the nutrient and 
chemoattractant helds, where the non-flux conditions are used at the chan¬ 
nel edges x=t) and L, i.e., {dS/dx)o = {dS/dx)i^ = 0 and {dN/dx)o = 
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Table 2: Parameter values used in simulations. 


Model parameters 


degradation rate of chemoattractant a 

5x10-3 [1/s] 

(0,2) 

production rate of chemoattractant b 

4x10^ [1/cell/s] 


consumption rate of nutrient c 

5 X 10-3^ [mL/cell/s] 

(1.0) 

diffusion coefficient of nutrient Dn 

8 X 10“® [cm^/s] 

(0.032) 

diffusion coefficient of chemoattractant Ds 

8 X 10“® [cm^/s] 

(0.032) 

mean tumble frequency ipo 

3.0 [1/s] 

(120) 

modulation of tumble freq. Xn/'^Po 

0.6 


modulation of tumble freq. Xs/i^o 

0.2 


coefficient in Eq. flTT]) cxi 

0.85 


coefficient in Eq. (ITT]) ct 2 

0.40 


stiffness of the response functions 5“^ 

8[s] 

(0.2) 

mean doubling time r 2 = In 2/r 

1.15 [h] 

(103.5) 

channel length L 

1.8 [cm] 

(18) 

Numerical parameters 


mesh interval Ax 

25 [/im] 

(0.025) 

time step size At 

0.2 [s] 

(0.005) 

total number of initial simulation particles 

56640 

(m;o=0.1) 


The values inside the parentheses indicate the values for the non-dimensional forms. 
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{dN/dx)j^ = 0 in Eqs. (12^ and fl25D . The boundary conditions for the sim¬ 
ulation particles are periodic in the y and z directions and specular reflection 
at the channel edges x=0 and L. The initial conditions for S and iV at t = 0 
are uniform and given as S{x) = 0 and N{x) = 1. The simulation particles 
are exponentially distributed according to the initial density prohle p{x), i.e., 

p{x) = aexp{—/3x), (35) 

where a and /3 satisfy the conditions p{x)dx = 1 and p{x)dx = 0.99 
with w = 2 [mm], 

4 . 1 . Macroscopic transport 

Figure |2] shows the traveling population wave of bacteria along the chan¬ 
nel. The population wave propagates with a constant velocity f^ave=0.16, 
i.e., Vwave=4.0 [/im/s], which is close to the experimental measurement of 
hwave=4.1 [/im/s]. The wave prohle changes only slightly after the initial 
transient period, t > 50; the density at the peak and rear side of the wave 
prohle only slightly increases as time progresses, whereas the front-side prohle 
of the wave is almost unchanged. 

The concentration prohles of chemical cues, i.e., the nutrient N and the 
chemoattractant S, are coupled to the collective motion of bacteria via the 
consumption of nutrients and the production of chemoattractants by bacte¬ 
ria. The motions of bacteria are also ahected by the concentration prohles 
of chemical cues. Figure [3] shows the concentration prohles of chemical cues 
and the spatial variation of the mean tumbling rate < T > along the relative 
coordinate x*, which is the position relative to the peak position that moves 
with a constant traveling speed f^ave- The local mean tumbling frequency is 
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Figure 2: The traveling population wave of bacteria along the channel, (a) Snap shots 
of the density profiles of the bacteria and (b) the superposition of the snap shots in (a) 
using the relative coordinate i*, which is the position relative to the peak position of the 
density profile that moves with a constant traveling speed t^ave- The traveling speed is 
calculated as V^ave=0.16, i.e., t4'ave=4.0 /rm/s. The dashed line in Fig. (a) shows the 
initial density profile of the bacteria. See also Supplemental Material for 3D motions of 
bacteria and time progress of bacterial density and nutrient concentration. 

obtained by taking an ensemble average of the tnmbling freqnencies of each 
bacterinm (see Eq. fl28|l i within the lattice site. In Fig. [3l every profile 
is also time-averaged over a time period t =50 to t =100. Incidentally, in 
Fig. [3l the mean tnmbling rates < > are not shown in the region x* > 0.7 

becanse the bacteria are not in this region. 

The nntrient concentration exponentially increases as the relative coordi¬ 
nate X* for X* < 0.2 and approaches nnity for x* > 0.5. The gradient takes a 
maximnm aronnd the position that coincides with the peak of the bacterial 
density, i.e., x* ~ 0, where the consnmption of nntrient by bacteria concen¬ 
trated. In contrast, the concentration of chemoattractant is not a monotonic 
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Figure 3: The spatial variations of the concentration of chemical cues, i.e., the nutrient 
N (solid line) and the chemoattractant S (dashed line), and the mean tumbling rate of 
the bacteria < ^ > (dotted lines) along the relative coordinate x* are shown. (See also 
the caption of Fig. H) The upper blue and lower orange dotted lines show the mean 
tumbling frequency of bacteria with negative and positive velocities, respectively, relative 
to the traveling speed e* (i.e., e* = — l^ave)- The left and right vertical axes show the 

chemical cues and the mean tumbling frequency, respectively. 

function; it moderately increases as x* at the rear side of the population 
wave, i.e., < 0 and takes a maximum around the peak of the density of 

bacteria, a;* ~ 0; thus, the concentration prohle of chemoattractant is similar 
to that for the bacterial density because the production of chemoattractant 
is proportional to the bacterial density. The chemoattractant concentration 
also exponentially decreases as x* due to a simple diffusion process in the 
region ahead of the population wave. 

The tumbling frequency of a bacterium depends on the material deriva¬ 
tives of chemical cues along the pathway of the bacterium. See Eq. fl28|) . 
Because the prohles of chemical cues do not change much along the relative 
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coordinate £*, the material derivative can be estimated as 


D 

'Dt 




(36) 


Thus, in Fig. |3l the modulation amplitudes of mean tumbling frequency 
for each of the positive and negative relative velocities are almost symmetric 
about the basal mean tumbling frequency T = 1. The modulation amplitudes 
are magnihed just behind the peak of the bacterial density, i.e., x* ~ 0, where 
the gradient of the nutrient takes a maximum and that of the chemoattrac¬ 
tant is also non-negative; the modulation amplitudes decrease as x* increases 
at the front side of the population wave a;* > 0, where both of the gradients 
of the chemical cues decrease. 


4-2. Microscopic dynamics 

The microscopic dynamics of the bacteria that form the traveling popu¬ 
lation wave are also investigated in terms of the probability density function 
(PDF) and the autocorrelation function (ACF) of the velocity of bacterium. 
Figure m shows the local PDFs p{ea) at different positions along the relative 
coordinate x*. The local PDFs are calculated by taking the time averages 
of the distribution of instantaneous velocities of the bacteria in each local 
lattice site for t=[50,100]. Note that the simulations are performed with a 
weight tCo=0.01, to obtain accurate prohles of the PDFs. The PDF of the 
longitudinal velocity, p{ex), rapidly increases around the traveling speed of 
the population wave, Cx ~ f^ave- The gradient around the traveling speed is 
larger in the region around the peak of the population wave, i.e., a;*=0 and 
±0.3, than in the rear and front regions, |£*| > 0.3. In contrast, the PDF of 
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the lateral velocity, p{ey), is symmetric and rather flat except at large veloc¬ 
ities, e-y ^ 1. The spatial variation of the PDF of the lateral velocity is also 
small. The steep gradient of the velocity distribution around the traveling 




Figure 4: The probability density functions of the velocity of a bacterium p{ea) at 
different positions along the relative coordinate x*. (a) The PDFs for the longitudinal 
velocity and (b) for the lateral velocity. The downward arrow shows the traveling speed 
of the population wave. 


speed in Fig. HJ^a) is related to the stiffness of the response function and the 
large spatial gradient of the nutrient around the peak. A detailed discussion 


is given in Appendix B 


The ACF of the deviation velocity, which is defined in Eq. fl38|l . of the 
bacteria that form the traveling population wave are shown in Fig. [5l The 
ACF G{f) is calculated from the trajectories of the velocity of each bacterium 
within a concentrated region, where the local density is not less than 10 % 
of the peak density at f=100, i.e.. 


G'(^) = - f)), {a = x,y,z), 


(37) 
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where 


= ea{i){i) - {ea{i){i)) . (38) 

Here, (•) represents the ensemble average of the test particles, and A(t) 
represents the time average of A{i) over i = [50,100]. The power spec¬ 
trum G{f) is calculated using the Fourier transform of G{t) as G{f) = 
Jq ° G(f) exp(—-i27r/f)df. Note that the inset in Fig. [5][a) magnihes the 
behavior of the ACF in a short time period scaled by the mean tumbling 
frequency -00 • 




Figure 5: The autocorrelation function G(r) (a) and power spectrum G{f) (b) of the 
deviation velocity, as defined in Eq. (I.SSI) . of bacteria within the traveling population 
wave. In the figures, t is the lag time, and / is the frequency, i.e., 1// is the period. The 
inset in (a) magnifies the behavior of the autocorrelation function in a short time period 
rescaled by the mean tumbling frequency ipo- 

The ACFs of the lateral deviation velocities, i.e., a=y and z, exponen¬ 
tially decrease on the short time scale, so that the autocorrelation of the lat¬ 
eral deviation velocities almost vanish after several tumbling events. These 
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features of the lateral deviation velocity demonstrate the simple Poisson pro¬ 
cess of the tumbling events of bacteria for movements in lateral directions. 

However, the ACF of the longitudinal deviation velocity exhibits consid¬ 
erably different behavior from that of the lateral deviation velocity. Notably, 
a negative autocorrelation arises for the longitudinal deviation velocity after 
a rapid decline on the short time scale due to the Poisson tumbling pro¬ 
cess. The amplitude of the negative correlation is small, but the negative 
correlation extends for a long time period (3> l/V’o)- The power spectrum 
of the longitudinal deviation velocity also shows a decline in the long time 
period. Thus, the negative autocorrelation function and its power spectrum 
characterize the random motions of the bacteria trapped within the popu¬ 
lation wave. The time period of the negative autocorrelation, fp, can be 
estimated as the characteristic time in which the traveling population wave 
passes through the width of density prohle Wp, estimated as Wp ~ Lq, i.e., 
"fp l/hwave- 

4 . 3 . Effect of variations in the parameters of the response function 

In the present kinetic model, the microscopic dynamics of each bacterium 
is coupled to the macroscopic transport of chemical cues via the response 
function dehned in Eq. ([7]), in which the sensitivity of the bacterium and 
the modulation amplitude of the tumbling frequency are characterized by 
the stiffness parameter and the modulation parameters xn and xs, re¬ 
spectively. The solutions obtained using the present kinetic model are signih- 
cantly affected by those parameters. As the stiffness parameter 5“^ increases, 
the prohle of the response function becomes similar to a step function; 

thus, the tumbling rate of the bacterium switches as soon as the sign of the 
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gradient of the chemical cue along the pathway changes. As the modula¬ 
tion parameter y increases, the difference in the mean tumbling frequencies 
of the bacteria for the negative and positive run velocities becomes larger; 
thus, the biased motion of the bacterium toward the migration direction is 
enhanced. In this subsection, the effects of changing the parameters of the 
response function are investigated. Monte Carlo simulations are performed 
for various values of the stiffness parameter, i.e., (5“^=0.125, 0.15, 0.2, 0.25, 
0.5, and 1.0, and the modulation parameter of nutrient, i.e., XAr=0.5, 0.55, 
0.6, 0.65, 0.7, and 0.8. 

Figures |6] and [7] show the traveling speed and density prohles of pop¬ 
ulation waves for different values of the stiffness parameters 5“^ and the 
modulation parameter of nutrient x^, respectively. The stiffness parameter 
of the response function 5“^ does not substantially affect the traveling speed 
but signihcantly affects the symmetry of the wave prohle. See Figs. Ela) 
and[71^a). As the stiffness increases, the width of the density prohle becomes 
thinner, and the peak of the density increases. The symmetry of the density 
prohle also changes. For small values of stihness, the tail created behind the 
population wave extends to a broader range. However, for the large values 
of stihness, i.e., 5“^=0.5 and 1.0, the density declines very steeply behind 
the peak, so the front side of the population wave becomes broader than the 
rear side. The traveling speed of the population wave Kvave shows a non¬ 
monotonic dependency of the stihness. The traveling speed is maximized at 
(5“^=0.2 in the present simulations. The mechanism for the maximum arising 
in the traveling speed with variations in the stihness parameter is not clearly 
known. However, it seems to be related to the relationship between the sen- 
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sitivity of bacterium to the nutrient and the concentration prohle of nutrient 
produced by the bacteria; as the stiffness parameter increases, the response 
of bacterium to the nutrient becomes more sensitive, but the width that the 
bacteria climb coherently along the gradient of nutrient become narrower. 

In contrast, the modulation parameter for nutrient in the response func¬ 
tion xtv does not substantially affect the prohle of the population wave. The 
symmetry of the wave prohle does not change, but the peak of the density 
decreases inversely to the length of the tail behind the population wave as the 
modulation parameter xn increases. However, the traveling speed is linearly 
related to the modulation parameter: l^ave oc xn- This linear dependency 
is also observed in the results obtained by the analytical formula of the trav¬ 
eling speed for the stepwise response function in the continuum limit, which 


is obtained in Ref. 


28| . The discrepancy between the results of the present 


Monte Carlo simulations and analytical formula arises because the present 


results are obtained for the hnite va 
stihness parameter 5“^ while, in Ref. 


nes of the Knudsen number and 


28| . the limiting case for those param¬ 


eters, i.e., —)■ 0 and 5~^ —)■ cxo, is considered. The convergence of the 

traveling speed obtained by the Monte Carlo method in that limiting case 
can be seen in Fig. [T31 

The prohles of the concentration of chemical cues, N and S', and the 
spatial variation of mean tumbling frequency < T >, which are shown in Fig. 
[31 are not substantially affected by changes in the modulation parameter yjv- 
The stiffness parameter 5“^ also does not substantially affect the prohles in 
the range of (5”^=0.125 to 0.25. However, for large stiffness parameter, i.e., 
5“^ =0.5 and 1.0, the prohles are signihcantly changed from the prohle shown 
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(a) (b) 

Figure 6 : The traveling speed of the population wave vs. the stiffness parameter 

(a) and the modulation parameter xtq (b) of the response function. The squares □ show 
the present results. Note that, for the present results, the modulation parameter is fixed 
as X 7 v= 0.6 in figure (a), the stiffness parameter is fixed as S~^=0.2 in figure (b), and the 
mean tumbling frequency is fixed as '0o=12O in both figures. The triangles A in figure 

(b) show those obtained by the analytical formula in Ref. for the stepwise response 
function (i.e., <5“^ —>■ oo) in the continuum limit where the Knudsen number vanishes. 







Figure 7: The superposition of the density profiles at i=75 along the relative coordinate 
X*” with variations in the stiffness parameter i5“^ (a) and modulation parameter of the 
nutrient xn (b). 

in Fig. |3l Figure E] shows the spatial variations of the chemical cues and the 
mean tumbling frequency along the relative coordinate x* for (5“^=1.0. As 
shown in Fig. [TJ the density of bacteria for (5“^=1.0 steeply declines behind 
the peak of the population wave, and the front side of the wave broadens. 
The concentration prohle of the nutrient N and the spatial variation of the 
mean tumbling frequency < T > for (5“^=1.0 are quite different from those 
in Fig. [3l In comparison with that in Fig. |3], the gradient of the nutrient N 
for (5“^=1.0 is rather flat behind the population wave, i.e., x* < —0.2, and 
is much steeper in the vicinity of the peak of the population wave, x* ~ 0. 
The gradient of the chemoattractant S is also larger behind the peak of the 
population wave, x* < 0. These prohles of chemical cues generate a local 
peak in the modulation of the mean tumbling frequency < T > behind the 
peak of the population wave. 
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Figure 8: The spatial variations of the concentration of the nutrient N (solid line) and the 
chemoattractant S (dashed line) and that of the mean tumbling rate of bacteria < 'I' > 
(dotted lines) along the relative coordinate x* for the stiffness parameter (5“^=1.0. See 
also the caption in Fig. |3] 

The effects of changing the stiffness and modnlation parameters of the 
response fnnction on the microscopic dynamics of bacterinm are shown in 
Figs. [9]and[T0l Fignre|9](a) shows the effect of changing the stiffness param¬ 
eter of the response fnnction on the PDF of the longitndinal velocity of 
bacterinm. As the stiffness parameter increases, the gradients of the PDFs 
in the vicinity of the traveling speed, ex ~ 0.15, become larger. The PDF 
for the large stiffness parameter, (5“^=0.5 and 1.0, shows an overshooting 
behavior after the steep gradient in the vicinity of the traveling speed. This 
overshoot profile is related to the asymmetric profile of the popnlation den¬ 
sity aronnd the peak ~ 0 (See Fig. [7]) and the stepwise response fnnction 
dne to the large stiffness parameter and steep gradient of the nntrient aronnd 
the peak. An intnitive explanation of the overshoot profile is also given in 
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Appendix B For a rigorous mathematical description, one can refer Ref. 


27|. 


In contrast, changing the modulation parameter xn does not affect the 
gradient of the PDF in the region around the traveling speed, but, at the 
left and right sides of the steep gradient region, the value of the PDF de¬ 
creases and increases, respectively, as the modulation parameter increases. 
(See Fig. [9](b).) The variation of the PDF at large longitudinal velocities due 
to changes in the modulation parameter is larger than that at small longitu¬ 
dinal velocities, so the mean velocity increases as the modulation parameter 
increases. 




Figure 9: The probability density functions of the longitudinal velocity of bacterium at 
the peak of the population wave i*=0 with variations in the stiffness parameter <5“^ (a) 
and modulation parameter xn (b). 


Figure [10] shows the ACF of the deviation velocity in the longitudinal 
direction of the bacteria that form the traveling population wave for various 
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values of the stiffness parameter of the response function. The ACF is only 
slightly affected by the modulation parameter but is signihcantly affected by 
the stiffness parameter. This fact indicates that the ACF, which represents 
the microscopic motions of bacteria in the traveling population wave, is highly 
related to the prohle of the population wave. (See also Fig. [71) As the 
stiffness parameter increases, the amplitude of the negative correlation, which 
arises just after an exponential decline due to the Poisson tumbling process, 
becomes larger. It seems that this behavior is related to the decline of the 
traveling speed at the large stiffness parameters shown in Fig. |6l 
0.002 


0.001 

G(f) 

0 


- 0.001 

Figure 10; The autocorrelation function G(f) of the longitudinal relative velocity of the 
bacteria that form the population wave for different values of the stiffness parameter 5~^ 
(a) and the modulation parameter (b). See also the caption in Fig. |SJ 

In summary of this subsection, the stiffness parameter 5“^ of the response 
function, which represents the sensitivity of a bacterium to chemical cues, 
signihcantly affects the wave prohle as well as the microscopic motions of bac¬ 
teria. The traveling speed is only slightly ahected by the stihness parameter. 
In contrast, the modulation parameter is linearly related to the traveling 
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speed but does not have a significant relationship with the wave profile or 
the symmetry of the velocity distribution function of the bacterium. 


4.4- Accuracy test 

The accuracy tests are conducted for various numerical parameter settings 
listed in Table |3l The STD corresponds to that listed in Table [2] and is used 
in the results shown in the former subsections. The T1~T3, M1-M3, and 
W1-W3 are used to investigate the effects of the variations in time-step size 
At, mesh interval Ax, and weight of a simulation particle Wq, respectively, on 
the convergence of density profile p{x) and traveling speed f^ave- The errors 
of the density profile and traveling speed are estimated as, respectively, 

Errp = i [ \p'{x) — p"{x)\dx, (39) 

L Jo 


and 


Err 


IE' 

I * W£ 


V" 


vw. 


V 

* wave 


( 40 ) 


where the prime and double prime “"’’marks indicate the coarser and 
finer numerical parameters, respectively. 

Figure fTTl shows the snapshots of the density profiles at t = 60 obtained 
for different parameter settings listed in Table |3l Table 0] shows the results 
of the error estimates Eqs. (I39|) and (ITOD . 

Except for Tl, the variation in time-step size At affects the density profile 
only in the vicinity of the peak but less affects the traveling speed. How¬ 
ever, only Tl gives a completely different solution. This is because for Tl, 
the product of the time-step size and tumbling frequency may 

exceed the unity. In the present Monte Carlo method, the tumbling of each 
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Table 3: The numerical parameter settings used for the accuracy tests. 



At 

Ax 

Wo 

STD 

0.005 

0.025 

0.1 

T1 

0.075 

- 


T2 

0.0025 

- 

- 

T3 

0.001 

- 

- 

Ml 

- 

0.1 

~ 

M2 

- 

0.05 

- 

M3 

- 

0.0125 

- 

W1 

- 


4.0 

W2 

~ 

- 

1.0 

W3 

~ 

- 

0.05 


The bar indicates the same value as for the STD. 
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(a) 



(b) 



(c) 


Figure 11: The snapshots of the density profile at t = 60 obtained for the different 
numerical parameters in Table [31 
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Table 4: Error estimates between different numerical parameter sets. 


<1 

Errp 

Err,> 

V wave 

STD-T3 

0.23 

3.8 X 10-^ 

T2—T3 

0.12 

2.3 X 10-^ 

(T1 gives 

a completely different profiles.) 

Ax 

Errp 

Err,> 

l^wave 

Ml-STD 

0.65 

3.0 X 10-2 

M2-STD 

0.24 

8.1 X 10-3 


(M3 induces a numerical divergence.) 


Wo 

Err^ 

Err,> 

wave 

W1-W3 

0.51 

3.8 X 10-3 

W2-W3 

0.21 

1.9 X 10-3 

STD-W3 

0.074 

3.4 X 10-^ 
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simulation particle, say the /tli particle, is stochastically determined with a 
probability ■ Thus, the time-step size At must be smaller than 

the inverse of the maxium tumbling frequency V'Max, which can be estimated 
as 

= V'o X (l + . (41) 

The mesh interval Ax affects both of the density prohle and traveling speed. 
The computation with M3 setting induces a numerical divergence because 
this numerical parameter setting does not satisfy the diffusion condition of 
the transport equations of chemical cues. At < 4^^ . Thus, it is important 
to note that there are two critical conditions for the time-step size At and 
mesh interval Ax to be satished, i.e.. 


At < -J—, 

(42) 

V’Max 

Ax‘^ 

(43) 

At < . . 

‘^Ds,n 


If the above conditions are satished, the numerical solution converges as 
setting the numerical parameters hner. 

The weight parameter wq (or the number of simulation particles) concerns 
the huctuation in the instantaneous density prohle, but the traveling speed 
is less ahected by the weight parameter. The amplitude of huctuation mea¬ 
sured by Eq. (I3^ decreases approximately in proportion to the square root 
of weight parameter Wq although the computational cost becomes larger 
inversely proportional to the weight parameter wq. 
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5. Comparison with asymptotic solution 


In this section, the results of the Monte Carlo simulations are compared 
with the macroscopic transport equation of the population density of bac¬ 
teria, which is obtained by the asymptotic analysis of the kinetic transport 
equation. Here, we only consider the case of the uniform scattering kernel, 
K = l/(47r), and nonproliferation, f = 0, for the kinetic transport equation. 

By introducing a new characteristic time dehned as = ^lJoLl/VQ, the 
kinetic transport equation Eq. (lT6l) for the uniform scattering kernel and 
nonproliferation is rescaled as 


dt dx a 


dvre 


^(e')/(t, X, e')dil{e') — 47r^(e)/(t, x, e] 


' all e' 


(44) 


where e corresponds to the Knudsen number, which is dehned as the ratio of 
the mean free path of bacterium Vq/V’o to the characteristic length Lq of the 
system and written as the inverse of the non-dimensional parameter of the 
mean tumbling frequency, i.e., e = (See also Sec. 12.21 ) Hereafter, the 
superscript “tilde”, e.g., t, represents the non-dimensional variable rescaled 
with the characteristic time 

Suppose that the modulations and xs are of amplitude e, respectively, 
i.e., Xn,s = ^4’n,s in Eq- flT8|) . and consider the asymptotic limit (or the 
continuum limit) e —?• 0 for Eq. fl44)) . Then, we obtain the following drift- 
diffusion equation for the population density of bacteria p(t, a:),|22| (See also 
Appendix) 


dp 

di 


= D. 


d^p d 


+ 


dx‘i dx. 


p{Ua[hgS] Ma[logiV]) 


(45) 
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where the flux is the functional of the gradients of chemical cues and 
defined as 

“«|r] = -^^£etatih (i|VF|{j d(. (46) 

The velocity distribution of bacteria becomes uniform in the continuum limit 
e —)■ 0, i.e., /(t, X, e) = p{t, x)/Att. Note that the reaction-diffusion equations 
of chemical cues Eqs. (na and ffTSjl are also rescaled by the characteristic 
time tg. The relations between the non-dimensional forms for the diffusion 
coefficients of chemical cues Ds^n-, the production rate of attractant a, the 
consumption rate of nutrient c, and time t scaled by tg (described by the 
superscript “tilde”) and to (described by the superscript “hat”), respectively, 
are written as 

Dn,s = Dn,s/^, a = a/e, c = c/e, t = et. (47) 

The Monte Carlo simulations are conducted with variations in the Knud- 
sen number e, i.e., e=0.02, 0.01, 0.005, 0.002, and 0.001, whereas the values of 
the other non-dimensional parameters are fixed as shown in Table O These 
values are set to coincide with those listed in Table [2] at '00 = 120. The drift- 
diffusion equation for the population density of bacteria p, Eq. flTSl) . coupled 
with Eqs. (1141) and f[T4D . is numerically calculated with the finite volume 
method for the parameters listed in Table [5l The flux is calculated at 
each time step with using the Simpson’s integral method according to the 
gradients of chemical cues obtained at the previous time step. To obtain 
accurate results for small Knudsen numbers in the Monte Carlo simulations, 
which are compared with the asymptotic solution in the continuum limit, the 
time-step size and particle number are set as At = 1 x 10“^ and Mg = 226560 
(or wq = 0.025), respectively, in the Monte Carlo simulations. 
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Table 5: The values of non-dimensional parameters for asymptotic equations. 


degradation rate of chemoattractant a 24.0 

production rate of chemoattractant b 1.0 

consumption rate of nutrient c 120.0 

diffusion coefficient of nutrient Dn 3.84 

diffusion coefficient of chemoattractant Ds 3.84 
modulation of tumble freq. (j)^ 72.0 

modulation of tumble freq. (j)s 24.0 

stiffness of the response function 0.2 


Figure O shows the snapshots of the population density of bacterial p{x) 
at t = 0.5 for various Knudsen numbers, i.e., e=0.02, 0.013, 0.01, 0.005, and 
0.001, obtained by the Monte Carlo simulations and the snapshot for the 
asymptotic limit obtained by the hnite volume calculation of Eqs. 045 p and 
fl4^ . coupled with Eqs. (IT^ and OT^ . It is seen that the results of the 
Monte Carlo simulations are asymptotically close to the asymptotic solution 
in the continuum limit as the Knudsen number decreases; the snapshot of 
the population density for e = 0.001 obtained by the Monte Carlo simulation 
almost coincides with the asymptotic solution in the continuum limit. How¬ 
ever, signihcant deviations from the asymptotic solution are also observed for 
moderately small Knudsen numbers, i.e., e >0.01. These deviations demon¬ 
strate that the asymptotic solution in the continuum limit is no more valid 
for that moderately-small Knudsen number. Figure [13] shows the conver¬ 
gence of the traveling speed Kvave scaled by e in the continuum limit e —)■ 0. 
It is seen that as increasing the stiffness parameter 5“^, the results obtained 
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Figure 12: Comparison of the population densities of bacteria for different Knudsen num¬ 
bers and for the asymptotic limit at a time t = 0.5. 


with Eq. (|T8ll converges to those obtained with the sign response function 
and the result of the sign response function converges to that obtained by 
the analytical formula in Ref. 22| in the continuum limit e —)■ 0. This also 
demonstrates that the present Monte Carlo method can accurately reproduce 


the analytical result of the traveling speed obtained in Ref. 


22 |. 


Figure [m shows the deviations of the probability density functions of lon¬ 
gitudinal velocity p{ex) from the uniform probability density po scaled by the 
Knudsen number e at the peak of population density x* = 0, which repre¬ 
sents the first order probability density function in the asymptotic expansion 
of e. It seems that the first order probability density asymptotically con¬ 
verges as the Knudsen number e decreases from 0.02 to 0.002, although the 
convergence of the hrst order probability density in the asymptotic limit is 
not very clearly seen because the fluctuation is notable at e = 0.001. This re- 
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Figure 13: The traveling speed of the population wave vs. the inverse of the Knudsen 
number for different stiffness parameters. For 5“^ —>■ oo, the hyperbolic tangent in Eq. 
(fT^ is replaced with the sign function X/|X|. The left arrow shows the traveling speed 


for 6 


oo in the continuum limit obtained by the analytical formula in Ref. 
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suit also indicates that the probability density p{ex) converges to the uniform 
probability density Pq in the continuum limit e —)■ 0. 

It should be noted that the Knudsen number depends not only on the 
biological properties of bacterium but also on the characteristic length Lq 
of the system. Thus, the Knudsen number becomes significant when one 
considers, for example, the cluster formation of bacteria in the micro devices 
and the collective migration of bacteria in the vicinity of a boundary. The 
deviations observed for moderately-small Knudsen numbers indicate that the 
present Monte Carlo method can take on a signihcance in the investigation 
on the micro scale systems. 

In the summary of this section, the asymptotic behaviors of the pop- 
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Figure 14: Comparison of the first-order probability density functions of longitudinal 
velocity for different Knudsen numbers at the peak of the population wave, po is the 
uniform distribution, i.e., po = 0.5, and represents the probability density in the continuum 
limit e —>■ 0. 

Illation density of bacteria and the velocity distribution of the bacterium 
are numerically demonstrated and the results obtained by the Monte Carlo 
simulations for small Knudsen numbers are compared to the asymptotic so¬ 
lution in the continuum limit. The comparison to the asymptotic solution 
validates the accuracy of the present Monte Carlo method in the asymp¬ 
totic behaviors for small Knudsen numbers. The distinctive deviations of the 
results obtained by Monte Carlo simulations from the asymptotic solution 
are observed in both the population density and the velocity distribution of 
bacteria for moderately small Knudsen numbers e >0.01. This result demon¬ 
strates the restriction of the continuum model and the utility of the present 
Monte Carlo simulation in the investigation on the micro-scale systems with 
moderately-small Knudsen numbers. 
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6. Summary 


In this paper, a Monte Carlo simnlation method for chemotactic bacteria 
is newly developed on the basis of the kinetic chemotaxis model, which was 
proposed in Ref. 28| . In this method, the Monte Carlo algorithm is employed 
to calculate the microscopic motions of bacteria that can be described using a 
model response function and a scattering kernel. The response function and 
scattering kernel depend on the local concentration of chemical cues, which 
is calculated using a hnite-volume method, so the microscopic motions of 
bacteria are coupled to the macroscopic transport of chemical cues via the 
response function and scattering kernel. 

The present simulation method can successfully reproduce the traveling 
population wave of chemotactic bacteria in a microchannel (Fig. [2]). The 
traveling speed and the wave prohle obtained by the simulation are close 
to those obtained in experiments. The microscopic dynamics of the bacte¬ 
ria that form the traveling population wave are investigated in terms of the 
probability density function (PDF) and the autocorrelation function (ACF) 
of the velocity of bacteria (Figs. IHandE]). The PDF of the lateral velocity 
of a bacterium is almost constant and uniform, except when the amplitude 
of the velocity is large. However, the PDF of the longitudinal velocity mono- 
tonically increases as the longitudinal velocity increases, so the PDF is biased 
toward the positive velocity regime. The PDF of the longitudinal velocity 
also spatially changes; the gradient of the PDF in the vicinity of the traveling 
speed of the population wave increases as the spatial position approaches the 
peak of the population wave, and the step-function-like prohle is obtained at 
the peak of the population wave. This stepwise prohle is related to the steep 
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spatial gradient of nutrients around the peak of the population wave. 

The ACF of the deviation velocity, which is calculated relative to the 
mean velocity, in the longitudinal direction shows a rapid exponential decline 
on the short time scale due to the Poisson tumbling process as well as the 
negative correlation that occurs after the Poisson tumbling process, which 
lasts a long time period. The negative autocorrelation function characterizes 
the random motions of the bacteria trapped within the population wave. 
The time period of the negative autocorrelation can be estimated as the 
characteristic time in which the traveling population wave passes through 
the width of density prohle. 

The response function, which determines the microscopic behavior of bac¬ 
teria, combines two important parameters, that is, the stiffness parameter, 
which represents the sensitivity of bacteria, and the modulation parameter, 
which describes the biased movements of bacteria toward chemical cues. The 
effect of changing the parameters of the response function on the macroscopic 
transport and the microscopic dynamics of bacteria is also investigated. The 
results showed that the stiffness parameter signihcantly affects both the wave 
prohle and the microscopic motions of bacteria. The traveling speed is only 
slightly affected by the stiffness parameter. In contrast, the modulation pa¬ 
rameter xn is linearly related to the traveling speed but does not signihcantly 
ahect the wave prohle or the microscopic motions of bacteria. 

The numerical accuracy of the method is tested with variation in the time- 
step size At, mesh-interval Ax, and weight parameter wq. It is found that 
the time-step size must be smaller than the minimum mean run duration, 
i.e.. At < Otherwise the Monte Carlo simulation produces a physically 
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incompatible result. The time-step size must satisfy the diffusion conditions 
of chemical cues, i.e., At < Ax^/{2Ds,n)- The weight parameter affects the 
fluctuation in the instantaneous density prohle but less affects the traveling 
speed. The fluctuation decreases as decreasing the weight parameter (or 
increasing the particle number). 

The drift-diffusion equation for the bacteria population density is derived 
by the asymptotic analysis of the kinetic chemotaxis equation in the contin¬ 
uum limit, where the Knudsen number e, which is defined by the ratio of 
the mean free path of bacterium to the characteristic length of the system, 
vanishes. In the present study, the results of Monte Carlo simulations for 
small Knudsen numbers are numerically compared with the asymptotic so¬ 
lution. It is demonstrated that the results of the Monte Carlo simulations 
for very small Knudsen numbers, i.e., e < 0.005, almost coincide with the 
asymptotic solution in the continuum limit but the deviation becomes no¬ 
table for e > 0.01. Thus, it is concluded that the continuum equation is 
valid for only a very small Knudsen number, e.g., e < 0.005, but for a hnite 
value of the Knudsen number, say e > 0.01, one needs to utilize the kinetic 
chemotaxis model. The comparison of the Monte Carlo simulation and the 
asymptotic solution also demonstrates the validity of the present method 
in the asymptotic behaviors of both the macroscopic transport of bacteria 
population density and the microscopic dynamics of the bacterium. 

It should be noted that the present kinetic chemotaxis model cannot treat 
the multi-body steric interactions and long-range hydrodynamic interactions 
between the bacteria, although those are supposed to be signihcant in a highly 
concentrated population. However, the kinetic chemotaxis model takes on 
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great significance in the investigation of the multiscale mechanism inherent 
in the chemotaxis problems in a moderate concentration where the steric and 
hydrodynamic effects can be ignored. 

In the present study, the Monte Carlo simulation is applied to a simple 
and fundamental problem of chemotactic bacteria in order to examine the 
validity of the present method. The method also succeeds in numerically 
demonstrating an asymptotic relationship between a continuum model and 
a kinetic model for the present problem. There remain two important chal¬ 
lenges to extending the present Monte Carlo method for applications to more 
complicated transport phenomena involved in practical engineering and bi¬ 
ological systems. First, the extension is of the geometry of the problem. 
The present Monte Carlo method treats the velocity space in the three- 
dimensional Cartesian coordinate, but the spatial geometry is restricted in 
the one-dimensional system. It is important to extend the method to the 
multi-dimensional space problem and demonstrate the colony dynamics in 
multi-dimensional spaces. Second, the extension is of the response func¬ 
tion to involve the memory effect of bacterium along its pathway. This 
extension allows us to investigate the memory effects of bacterium in the 


colony dynamics. 
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4Q|, 


41 


42| The particle-based algorithm employed 


m 


the present Monte Carlo method is thought to be useful in incorporating the 
memories of each bacterium into the present Monte Carlo method. These 
extensions are important future works. 
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Appendix A. Asymptotic analysis 


In this appendix, the macroscopic transport equation of the population 
density of bacteria is derived by the asymptotic analysis of the kinetic chemo- 
taxis model Eq. fl4T|) according to Ref. [2^. The drift diffusion equation 


of the population density of bacteria is obtained in the continuum limit, 
where the Knudsen number e, which is defined by the ratio of the mean free 
path of bacterium Vq/V'o to the characteristic length Lq of the system, i.e., 
e = Vq/'iPqLq, vanishes. The basic procedure of the asymptotic analysis was 
put forward in Ref. 22|, but the functional form of the response function of 
bacteria used in the present study is only different from that in the reference. 

We seek the asymptotic solution of the bacterial density f(t, x, e) in a 
expansion of e, i.e., / = /o + e/i + e /2 + • • • ■ By substituting the expansion 
form into Eq. fl44|) and arranging the terms by the order in e, we obtain the 
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following expansion, 


. dfo dfi 

dXa I di ^dXa 




fo{e')dVt{e') -47r/o(e) 


' all e' 


(^i(e')/o(e') + fi{e'))dQ{e') - 47r/i(e) + 47r4^i(e)/o(e) 


all e' 

2^ 


+ Oie^) } , 


(A.l) 


where ^i(e) is obtained by the expansion of Eq. 03 and is written as 


<i-,(e) = :| tanh (^-^) + ^ tanh f 


(A. 2) 


Here, we suppose that the modulation parameters xs and xn are of ampli¬ 
tude e, respectively, i.e., xs,n = ^4’s,n- 

For the limit e —)■ 0, we obtain an uniform distribution as to the bacterial 
velocity, i.e., 

fo{t,x,e) = — -, A.3 

47r 

where po{i,x) is the population density of bacteria in the asymptotic limit 
and is obtained by po = / /o(e)dr2(e). Hereafter, unless otherwise stated, 
the integral as to the vector e is performed over all directions of e. By 
substituting Eq. flA.Sp into Eq. flA.ip and integrating Eq. flA.ip as to e, we 
obtain the following equation for po, 

dpo 


di 


+ V ■ = 0, 


(A.4) 


where the flux ji is dehned as ji = f efi(e)dQ(e). By integrating Eq. flA.ip 
multiplied by e as to e, we obtain the equation for the flux ji in the first 
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order of e, 


dvr 


Cr^epdfl^e) ) + -^po j ea'^i{e)dQ.{e). (A.5) 


The integral in the hrst term of the right hand side of Eq. flA.5l) is calculated 
as / eaepdVL^e) = (47r/3)5a/3. The integrals of the hrst and second terms of 
Eq. flA.2p , which are involved in the second term of the right hand side of 
Eq. dQ, are calculated by using the polar coordinates whose pole directions 
are set to the directions of gradient vector V log S and V log N, respectively. 
Thus, we obtain the equation for the hux ji, 


ji = -DpVpo + po 


<Ps viog^ 

[4 |Vlog^|i_i 
V log AT 

4 iVlogiVli., 


^ tanh 
.^tanh 


(l 


VlogSIf df 


(l 


IVlogiVIJ 


(A.6) 


where Dp = 1/3. By substituting Eq. flA.bp into Eq. flA.dp . we obtain the 
drift diffusion equation for the population density in the asymptotic limit. 


Appendix B. Overshoot profile of the velocity distribution 


In this appendix, we discuss the overshoot prohle of the velocity distri¬ 
bution which appears at the peak of the population wave for a large stiffness 
parameter 5“^. See Fig. [9l For the simplicity, we consider the uniform scat¬ 
tering kernel K = l/dvr and nonproliferation f=0 because these factors less 
affects the results. In the one-dimensional relative coordinate system as in¬ 
troduced in Eq. fl36|) . the density of the bacteria with a velocity e^,, /(x*, e^) 
can be written as 


G _ ( - Ewave ^ df 

Kex) V A(e,) J dx* 


(B.l) 
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where G is the generation rate due to the tumbling and is written as 


G 


1 

2 



K^x)f{x* ,ex)dea,, 


(B.2) 


and the tumbling rate X{ex) is written as 


A(e^) = '4’o 


1 Xn . 1 ( <9 log iV \ Xs , u ( ~ ^ S 

1 — — tcirin I ---- I — — tciriii I -:::— 


dx* 


■ 


dx* 


(B.3) 

Here, we assume the velocity space for e is axially symmetric and d/dy = 
d/dz = 0. Note that the generation rate G is independent of and A(ea;)“^ 
represents the mean run duration of the bacteria with a velocity ix- The 
mean run distance of the bacteria with a velocity Cx in the relative coordinate 
system can be written as 


nex) = 


lAvave 


A(e,^ 


The density of the bacteria with a velocity Cx at the peak of the wave x 
can be approximately written as, for a small I*, 


(B.4) 
* = 0 


/(o, Cx) = 


G 


1 r 


A (cj 


+ 2 


f(-r,ex)-f{t,ex) 


(B.5) 


Here, the hrst term of the right-hand side of Eq. fIB.Sp . which we call the 
generation term, represents the generation of a new velocity Cx due to the 
tumbling in a mean run duration A(ea;)“^ and the second term of the right- 
hand side of Eq. flB.Sp . which we call the flux term, represents the net flux 
of the density of the bacteria migrating without tumbling. 

The tumbling rate at the peak of the wave A(ea;) steeply decreases in the 
vicinity of Cx = f^ave because the spatial gradient of nutrient concentration 
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is very large at the peak of the wave. Thus, the velocity distribution /(O, e^) 
steeply increases in the vicinity of = f^ave due to the contribution of the 
generation term in Eq. fIB.SD while the contribution of the flux term is neg¬ 
ligible because the mean run distance vanishes at = Kvave- However, the 
flux term becomes signihcant as the velocity separates from the traveling 
speed V^ave- Especially, for a large stiffness parameter <5“^, the tumbling rate 
\{ex) is a stepwise function of Cx — Kvave, so that, except only the vicinity of 
= Kvave, the mean run distance l*{ex) is proportional to e^, — Kvave while 
the generation term does not vary signihcantly. In addition, the density pro- 
hle around the peak of the wave becomes highly asymmetric as the stiffness 
parameter 5“^ is large; the rear side of the wave steeply declines while the 
front side of the wave broadens. Thus, for e^, — Kvave > 0, the density at the 
front side of the peak /(/*, e^,) is larger than that at the rear side of the peak 
f{—l*,ex), and the flux term of Eq. fIB.Sp negatively grows as ix — Krave- 
This gives the overshoot prohle of the velocity distribution at the peak of the 
wave for a large stiffness parameter. 
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